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Planets in the habitable zone of lower-mass stars are often assumed to be in a 
state of tidally synchronized rotation, which would considerably affect their pu¬ 
tative habitability. Although thermal tides cause Venus to rotate retrogradely, 
simple scaling arguments tend to attribute this peculiarity to the massive Venu¬ 
sian atmosphere. Using a global climate model, we show that even a relatively 
thin atmosphere can drive terrestrial planets’ rotation away from synchronic- 
ity. We derive a more realistic atmospheric tide model that predicts four asyn¬ 
chronous equilibrium spin states, two being stable, when the amplitude of the 
thermal tide exceeds a threshold that is met for habitable Earth-like planets 
with a 1 bar atmosphere around stars more massive than ~ 0.5 — 0.1 Mq. Thus, 
many recently discovered terrestrial planets could exhibit asynchronous spin- 
orbit rotation, even with a thin atmosphere. 


As we all experience in our everyday life, atmospheric temperatures oscillate following the diurnal 
insolation cycle. This, in turn, creates periodic large-scale mass redistribution inside the atmosphere, 
the so-called thermal atmospheric tides. But as we all have experienced too, the hottest moment of the 
day is actually not when the Sun is directly overhead, but a few hours later. This is due to the thermal 
inertia of the ground and atmosphere that creates a delay between the solar heating and thermal response 
(driving mass redistribution), causing the whole atmospheric response to lag behind the Sun (i). 

Because of this asymmetry in the atmospheric mass redistribution with respect to the sub-solar 
point, the gravitational pull exerted by the Sun on the atmosphere has a non-zero net torque that tends 
to accelerate or decelerate its rotation, depending on the direction of the solar motion (2, 3). Because 
the atmosphere and the surface are usually well coupled by friction in the atmospheric boundary layer, 
the angular momentum transferred from the orbit to the atmosphere is then transferred to the bulk of the 
planet, modifying its spin (4). 
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On Earth, this effect is negligible because we are too far away from the Sun, but the atmospheric 
torque due to thermal tides can be very powerful, as seen on Venus. Indeed, although tidal friction 
inside the planet is continuously trying to spin it down to a state of synchronous rotation, thermal tides 
are strong enough to drive the planet out of synchronicity and to force the slow retrograde rotation that 
we see today (2-6). Very simple scaling arguments predict that the amplitude of the thermal tide is pro¬ 
portional to the ratio of the atmospheric mean surface pressure over its scale height (1). Everything else 
being equal, one would thus expect the thermal tide to be ~ 50 times weaker if Venus had a less massive, 
cooler Earth-like atmosphere. Does this scaling really hold and how massive does an atmosphere need 
to be to affect the planetary rotation, this has not been completely worked out yet. 

These questions are of utmost importance as we now find many terrestrial planets in a situation 
similar to Venus. Because the habitable zone is closer around lower-mass stars, planets in this region 
are often expected to be tidally synchronized with the orbit {7-14), which seems to create additional 
difficulties for their habitability. In particular, the permanent night side can be an efficient cold-trap for 
water {13, 14), strongly destabilize the carbonate-silicate cycle {10, 11), and even cause atmospheric 
collapse in extreme cases (9, 15). 




Eigure 1: Frequency dependence of the torque. 

Sin of the lag angle (sin26a; A), amplitude of the 
pressure bulge (|^a|; B), and torque (Im(^a); C) 
as a function of the forcing frequency {{(0 — n)/n) 
computed from the numerical atmospheric model 
(gray data points) and given by our analytical 
model (black curve). Results are shown for two 
pressures (1 and 10 bar) in the 1366 W/m^^ case 
with an orbital period of 225 d. The error bars 
show the internal variability (±la). The Venu¬ 
sian tide amplitude is shown in B. Despite its 
simplicity, the analytical model fairly captures 
the frequency dependence of the thermal tide 
response. 


Here, we investigate whether or not thermal tides can drive terrestrial planets in the habitable zone, 
with a relatively thin atmosphere, out of synchronous rotation. Previous studies on Venus (2-6), and for 
exoplanets {16,17), have shown that this reduces to the search for equilibrium rotation states for which 
the bodily torque (Eg) and the atmospheric torque (Ta) cancel each other and provide a restoring force 
against deviations from this equilibrium. In the circular case with zero obliquity, these torques are given 
by (5-6) 

3 3 

Ta = — -.fifa6a (2ft) — In) and Tg = — -Kgbg{l (0 — In), 
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where A'a = (3M*/?p/(5pa^), Kg = GM^R^^/a^, G is the universal gravitational constant, M* is the 
stellar mass, Rp, co, and p the planet’s radius, rotation rate, and mean density, and a and n its orbital 
semi-major axis and mean motion. bg{a) characterizes the frequency-(a-)dependent response of the 
body of the planet (i.e. its rheology), and ba{o) that of the atmosphere. Note that, unlike trapping in 
asynchronous spin-orbit resonances (like Mercury; 18), thermal tides do not need any eccentricity to 
drive a planet out of synchronous rotation. While some general conclusions can be reached without a 
precise knowledge of these responses (5, 6, 16, 17), the biggest limitation in predicting the properties of 
equilibrium spin states lies in the uncertainties on the shape and amplitude of Z>a(<7)- In particular, the 
relation between the mass of the atmosphere and the strength of the thermal tide, the key ingredient to 
quantify whether asynchronicity is ubiquitous, remains unknown. 

To tackle this issue, we use a generic global climate model {14, 19, 20), commonly used to model 
planets in the habitable zone of low-mass stars, to empirically quantify the amplitude of torque induced 
by thermal tides for planets with various atmospheric masses (characterized by the surface pressure, ps)> 
compositions, and incoming stellar fluxes {F). Once ps and F are chosen, we run the atmospheric model 
for several diurnal frequencies, a = (0 — n. Because of the thermal forcing, a surface pressure pattern 
lagging behind the substellar point forms (see Fig. SI). Once mean thermal equilibrium is reached, we 
compute the complex amplitude of the quadrupolar thermal tide {q^', 21), as shown in Fig. 1, and the 

value of the torque is given by equation (1) where ba{2(0 — 2n) = — ^J~^lm{qa{co — n)). 

To test our framework, we applied our model to the Earth and Venus, for which results meet existing 
constraints (see Fig. S2; 21). However, rather counterintuitively, for the same forcing frequency, the 
amplitude of the thermal tide in a 1 bar Earth-like atmosphere is almost an order of magnitude stronger 
than on Venus (Eig. 1 .B). This difference is the result of the sunlight being almost completely scattered 
or absorbed before it reaches Venus’ surface {21). 

As discussed earlier, the lag and amplitude of thermal tides are closely related to the thermal inertia 
of the system. In fact, a careful analysis of the result of the climate model shows that the frequency 
dependence of the atmospheric response shown in Eig. 1 is fairly analogous to the thermal response of a 
radiating slab with a finite thermal inertia that is periodically heated, so that, to a very good approxima¬ 
tion, we can write 


^a(Ct) 


9o 

l+ia/(Oo' 


( 2 ) 


where — 1, ft)o is the inverse of the timescale needed for the system to reach thermal equilibrium, 
and qo is the amplitude of the quadrupole term of the pressure field at zero frequency {21). In theory, 
COq can be estimated if the heat capacity of the system is known {21), but in practice, both coq and <70 are 
computed from the numerical model for a given atmosphere and are shown in Table 1 for limit cases. In 
the circular case with zero obliquity, the torque thus writes. 


Ta{o)-n,ps,F) 


-Kaqo{ps,F) 


{(0-n)/(Oo{p^,F) 

\ + [{co-n)/cOQ{p^,F)f 


(3) 


Interestingly, although governed by different physics, the atmospheric torque follows the same law as 
the body torque for a viscoelastic sphere (with the opposite sign). 

A first qualitative difference with previous works {3-5, 16, 17), is that we derive a very different 
functional form for the atmospheric torque. In particular, the function /(cj) = ba{2o)/bg{2o) (5) is not 
monotonic around potential equilibria when a realistic rheology is used. As seen in Eig. 2, for a constant- 
Q or an Andrade rheology, this results in the possible existence of up to five equilibria in the circular 
case, two of them being unstable {21). The diversity of equilibria might be even richer in eccentric 
systems where these numbers could change {16,17). Interestingly, the synchronous spin state is stable. 
Knowing that Venus, despite such a rheology, did not end up synchronized tells us that a planet can avoid 
being trapped in such a stable synchronous state and constrains the history of the Venusian atmosphere 
( 21 ). 
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Table 1: Numerical values of the amplitude of the atmospheric quadrupole (q'o) and intrinsic thermal frequency 
of the atmosphere (cOq) derived from the global climate model. 

Qf Model output 


simulations 

F {W.var^) 

Ps (bar) 

qo (Pa) 

®o(s ^) 

ln/cooid) 

Venus 

2610 

92 

201 

3.77x10-^ 

193 

Inner 

1366 

1 

1180 

2.30x10-® 

32 

habitable zone 


10 

4050 

1.46x10-® 

50 

Outer 

450 

1 

890 

1.18x10-® 

62 

habitable zone (N 2 ) 


10 

2960 

7.17x10-^ 

101 

Outer 






habitable zone (CO 2 ) 

450 

10 

2590 

9.7x10-^ 

70 
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Figure 2: Equilibrium spin states of the planet. Atmospheric (dashed), gravitational (dotted) and 
total (solid) torque as a function of spin rate for two tidal models (Andrade: A, B, C; Constant-Q: D, 
E, E). Arrows show the sense of spin evolution. A, D: weak atmospheric torque, only one equilibrium, 
synchronous spin state exists (blue disk). B, E: bifurcation point (a = a^). C, E: the atmospheric torque 
is strong enough to generate four asynchronous, equilibrium spin states, two being unstable (red circles) 
and two being stable (blue disks; one is retrograde in the case shown). The synchronous spin state 
remains stable. 
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In addition, the number and location of equilibria undergoes a bifurcation as asynchronous spin 
states exist only when the amplitude of the thermal tide reaches a threshold. Thus, our results reveal the 
existence a critical distance (Oc) beyond which the planet can be asynchronous which, using a constant-2 
rheology, reads 


ac 



GM^pRl^ 

qoQ 


1/3 


(4) 


where ^2 is the Love number and Q the tidal quality factor (2i). Both Uc (Fig. 3) and the equilibrium 
asynchronicity (\(0 — n\ = COq ; Fig. S3) can be computed for various cases 

using Table 1 . The corollary is that, even without any triaxiality, planets on circular orbits for which 
atmospheric tides are too weak should be in exact spin-orbit resonance. 

More importantly, our results provide a robust framework for the quantitative assessment of the 
efficiency of thermal tides for different atmospheric masses without having to rely on scaling argu¬ 
ments calibrated on Venus. This is crucial as Venus thermal tides turn out to be relatively weak (see 
Fig. 1 .B). As can be seen in Fig. 3, Earth-like planets with a 1 bar atmosphere are expected to have a 
non-synchronous rotation if they are in the habitable zone of stars more massive than ~0.5 to 0.7 Mq 
(depending on their location in the habitable zone). This lower limit becomes < 03Mq for a 10 bar 
atmosphere. Interestingly, these limits are much less restrictive than the one obtained from our Venus 
model (purple line in Fig. 3). This realization required full atmospheric modeling. 
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Figure 3: Spin state of planets in the hab¬ 
itable zone. Each solid black line marks 
the critical orbital distance {a^, Eq. (4)) 
separating synchronous (on the left, red ar¬ 
rows) from asynchronous planets (on the 
right, blue arrows) for Ps = 1 and 10 bar 
(the extrapolation outside the habitable 
zone is shown with dotted lines). Objects 
in the gray area are not spun down by tides. 
The habitable zone is depicted by the blue 
region {14, 24). Gray dots are detected and 
candidate exoplanets. The error bar illus¬ 
trates how limits would shift when varying 
the dissipation inside the planet {Q ~ 100; 
21) within an order of magnitude. 


Atmosphere as thin as 1 bar being a reasonable expectation value given existing models and Solar 
System examples, our results demonstrate that asynchronism mediated by thermal tides should affect 
an important fraction of planets in the habitable zone of lower mass stars. This is especially true in the 
outer habitable zone where planets are expected to build massive CO 2 atmospheres (7). 

This has many implications. On one hand, all the threats to habitability created by the presence of 
a permanent cold, night side may not be as severe as usually thought. On the other hand, the habitable 
zone has been recently shown to be more extended for synchronous planets (72). Eor these objects, if 
the atmosphere is thick enough, the non-synchronous rotation that will ensue will thus come to limit the 
extent of the habitable zone around lower-mass stars. 

The thermal tide mechanism presented here does not only affect habitable planets, so that many other 
terrestrial bodies with significant atmospheres could potentially have asynchronous rotations, depending 
on their orbital location (see Eig. 3). With that in mind, observational methods able to constrain the 
rotation rate of exoplanets (22, 23) become more valuable, and could even be used to constrain the 
atmospheres of the latter. 
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Material and Methods 


We performed simulations of the thermal response of the atmosphere of terrestrial planets in the habit¬ 
able zone for a wide range of forcing frequencies (solar day duration) using the LMDZ generic global 
climate model (Sect. 1). In each of these simulations, the response is characterized in terms of the com¬ 
plex amplitude of the mass quadrupole {qa{(0 — n,ps,F)-, (Sect. 2)), which allows us to compute the 
atmospheric torque exerted by the star on the atmosphere. This model is tested against constraints from 
the Earth and Venus (Sect. 3). 

We then show that, by analogy with the simple problem of a thermally forced, radiating slab, we are 
able to model analytically the frequency dependence of the atmospheric torque (Sect. 4). This allows us 
to provide an analytical formula that uses only two free parameters {qo{ps,F) and coo{ps,F)) which can 
be calibrated from the simulations and are given in Table 1 for models bracketing the behavior expected 
in the habitable zone. 

1 Numerical atmospheric simulations 

To model thermal tides, we need to know the redistribution of mass in the atmosphere. As we will 
see, because of hydrostatic balance of the atmosphere, one only needs to know the spatial and temporal 
distribution of surface pressure. To model the evolution of the atmosphere, we use the Generic version 
of the LMDZ global climate model that has been specifically designed for exoplanet studies {14, 19, 
20). The details of the version used here can be found in ref. 14. 

For simplicity, we only model zero obliquity planets without topography, with a uniform surface 
albedo and thermal inertia (0.2 and 2000 S.I. respectively), on circular orbits. The atmosphere is assumed 
to be dry and composed of an Earth-like mixture of N 2 with 376ppmv of CO 2 , and radiative transfer is 
computed consistently using the correlated-k method. 

Despite the fact that the orbital distance of the habitable zone changes drastically for different stars, 
the flux defining its boundaries only changes by ~ 25% across the main sequence. Therefore, to re¬ 
duce the computational cost, we ran suites of simulations for two insolations: i) an Earth-like case 
{F = 1366 Wm^^) which is representative of the inner edge, and ii) a F = 450 W.m^^ case which is rep¬ 
resentative of the outer edge. These two sets of simulations should encompass the behavior expected 
for planets in the habitable zone. As will be shown hereafter, however, the effect of the incoming stellar 
flux on the distance at which asynchronicity sets in is rather small compared to other effect such as 
the effect of mean surface pressure. For each insolation, several mean surface pressures (mainly 1 and 
10 bar) were used, and for each of these pressures, many simulations were done to cover the relevant 
range in forcing frequency m — n (or equivalently the solar day duration, l%j{(a — n)). This represents 
almost a hundred simulations. Snapshots of the pressure variations created by the diurnal tide in the 
{1366 W.m^, 10bar} case for four different frequencies are shown in Fig. SI. 

As discussed in the main text, thanks to the carbonate-silicate cycle, CO 2 is thought to accumulate 
in the atmosphere near the outer edge to keep the surface unfrozen (7). Our F = 450 W.m^^ cases with 
a N 2 dominated atmosphere are indeed too cold for liquid water. We thus ran an additional suite of 
models for a 10 bar, pure CO 2 atmosphere. Results show that changing from a N 2 to a CO 2 atmosphere 
reduces slightly the amplitude of the tide because of the increased albedo and higher absorption of CO 2 
(see Table 1). This effect is however less important than the effect of the reduced flux. 

As described hereafter, the amplitude and lag of the thermal tide is computed from each simulation 
(Sect. 2). Then, for a given set of parameters {FjPs}; the frequency dependence of the thermal tide is 
fitted following the method described in Sect. 4. Final results are given in Table 1. 

2 Torque due to thermal tides 

Because of the non-spatially homogenous insolation that it receives, the atmosphere is continuously 
driven out of dynamical equilibrium and undergoes large-scale motions. This results in a non-spherically 
symmetric redistribution of its mass which thus can be gravitationally torqued by the star. To compute 
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this torque, we first express the gravitational potential created by the atmosphere at a point r = (r, 0,0) 
in space 


U,{r) = -G 


/ 


dm' 

r — r' 



PW) 

r — r' 


r'^ sin o'd6'd^'dr', 


where G is the universal gravitational constant, p(r') is the density at the location r' = +z, O', 0'), 

and integration is carried out over the atmosphere. Assuming that the atmosphere is thin (r' ~ Rp) and 
close to hydrostatic equilibrium, we have 


dp 

dz 


-gp ^ 



g 


g 


where p is the pressure, Ps{0', <p') its value at a point of the surface, and g the surface gravity. 
Expanding |r — r'| in terms of Legendre polynomials (P/) and integrating over / we get 


f/a(r) 



J P/(r •r')ps(0^0O sinO'dO'd^'. 


Using the addition theorem, we expand P/ (r • r') in spherical harmonics, which yields 


^ Att / J? \ r 

Ua(r) =-/E/-(e',f)p,(0',f)sin0'de'df 

g 2/ + 1 V r y J 


GRr 


I I 


471 /P 


g ;n^;2Z + lVr/ 


\ ^+1 


^ PTYne,^), 


where p'" = f F™* ps sin O'dO'd^', are the complex moments of the pressure field. 

The gravifafional forque applied by fhe sfar on fhe afmosphere (fhe opposite of fhe torque by fhe 
afmosphere on fhe sfar) is fhus given by Ta = M*r x VUa- For fhe sake of concreteness, we confine 
ourselves fo fhe zero obliquify case where fhe forque is perpendicular fo fhe planefary equafor and is 


Pa = M*r 



GM^Rp 

g 


CO 


I 


471 

2Z+1 



/ 


£ imuryne. 

m——l 


71 

2 ’ 


<!>*), 


where fhe dipole, Z = 1 ferm has been discarded as if is only a change in fhe center of mass (6). The 
colafifude and longifude of fhe subsfellar poinf (6* and (p^,) are measured in fhe coordinafe sysfem defined 
by fhe spin vecfor and an arbifrary meridian of origin (fhe choice of fhe origin is inconsequential as long 
as fhe pressure momenfs are compufed in fhe same coordinafe sysfem). Because, in our simulafions, fhe 
amplifude of fhe p'” decreases af higher I and Pp/r <C 1, we can keep fhe 1 = 2 terms only. Furfhermore, 
considering fhaf (0* = |, (/)*)= 0 and p^^Y^^ = (— 1 )^^^P 2 *F 2 ^*, we have 



where Im(^a) = Im{p 2 e^"^*) = jp^l sin(2(/)* — 2^a)^ where 2(/)a is minus fhe argumenf of fhe complex 

number p^. We also infroduce 5a = (/>* — (pa- Defining 5a(2a) = Iin(^a(<7)). we recover fhe 

equation in fhe fexf. qa can be seen as fhe complex amplifude of fhe mass quadrupole in fhe reference 
frame where fhe subsfellar poinf is faken as fhe origin of longifudes. In a fypical sifuafion wifh no 
eccenfricify, bofh \p\\ and 5a show only small flucfuafions around a consfanf (see Fig. 1). For higher 
rofafion rates, baroclinic waves sfarl to appear which increases fhe internal variabilify. 
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For illustrative purposes, we show the semi-diurnal component of the surface pressure field, i.e. 
p ^2 + pressure maximum lagging behind the star does so with 

a lag of more than 90° (in other words, the pressure maximum closest to the star is leading the motion 
of the substellar point) which is why atmospheric tides have an effect opposite to body tides. For these 
slow rotations, the diurnal component is still dominating the pressure distribution at the surface. This 
is the primary reason for the strong day/night side circulation seen in many 3D atmospheric models of 
tidally locked exoplanets (9, 12, 14), but the amplitude of this component becomes smaller at Earth-like 
rotation rates (see below). 

3 Application to Earth, Venus, and more massive planets 

Our most extensive source of data on thermal tides coming from the Earth, we have validated our nu¬ 
merical model by running a simulation for this specific case (accounting for the actual topography and 
albedo variations of the planet; see ref. 20). To compare with meteorological data, we compute the am¬ 
plitude and phase of both the diurnal and semi-diurnal component of the tidal pressure wave at each 
surface grid point. This is done by performing a Eourier analysis of the pressure time series output by 
the simulation (i). In this case, the spatial distribution of the semi-diurnal phase and amplitude predicted 
by our atmospheric model is fairly close to the data (see Eig. S2 to be compared with Eig. 2S.3 of ref. 1). 
Our semi-diurnal solar tide has a ~2mb amplitude and « 150° phase at the equator which is in good 
agreement with the meteorological observations (~1.2mb and « 150— 160°; 1) We also recover the fact 
that, for a fast rotating planet like the Earth, the diurnal component of the tides has a small amplitude 
compared to the semi-diurnal one 

Note, however, that for the Earth, both gravitational and atmospheric torques (created by the Sun) 
are too weak to bring the Earth to an equilibrium spin state over its lifetime. This should be true for 
planets in the habitable zone of stars above ~ 0.9Mq. These objects are thus not expected to be found 
in an equilibrium spin state (either synchronous or not). We also tested the model on Mars where the 
amplitude of the pressure variations at the surface were found to be in agreement with surface data (25). 
We do not show these results here because the Martian atmosphere is too thin (~ 6 mb) to significantly 
affect the spin of the planet. 

Because the amplitude of the thermal tide is greater in a 10 bar atmosphere than in a Ibar atmo¬ 
sphere, one might expect that the tide would be stronger in a Venus-like atmosphere, because it is even 
more massive. There is, however, no contradiction with our theory as we find that the amplitude of ther¬ 
mal tides does not increase continuously with surface pressure. Tests show that above a surface pressure 
of a few tens of bars (depending on the atmospheric composition), the amplitude of thermal tides starts 
decreasing again. This is caused by stellar energy being progressively absorbed higher in the atmosphere 
and the tide being damped by a growing quiescent layer beneath it. Eor very massive atmospheres, the 
system should slowly transition toward a regime where the solid surface has a negligible effect, although 
the thermal tide itself might still impact the spin of the planet (26). The exact mass of atmosphere at 
which this transition occurs needs to be investigated in more detail. 

Venus, in addition to its thick atmosphere, exhibits very reflective and absorbing clouds that strongly 
decrease the amount of light absorbed in the lower atmosphere, and thus the forcing. Indeed, only a few 
watts per square meters reach the surface. To take that into account, we use outputs of the Venus version 
of the EMDZ model (27) and compute the amplitude of thermal tides as described above. This yields 
an amplitude |^a| ~l-1.6mb for the semi-diurnal solar tide (as shown in Eig. l.B), and a torque that is 
roughly equal to the gravitational torque expected for Venus (4) for a quality factor of 2 ~ 100, which 
is in the range of what can be expected for a rocky mantle without oceans. This shows that, even in 
this well constrained case, our framework yields sensible results. To facilitate comparisons, this value 
of the quality factor for bodily tides is the fiducial value that we use when needed, but we also show the 
sensitivity of our results on this parameter that may vary from a planet to another (especially in Eig. 3). 
The values of qo and coo needed to reproduce are given in Table 1. 
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4 A thermal inertia model for the atmospheric bulge 


In first approximation, the atmospheric bulge created by the periodic insolation pattern is not aligned 
with it only because it takes some time for the atmosphere and surface to increase their temperature 
when heated. Thermal inertia of the atmosphere-surface system should thus be an important parameter 
in determining both the amplitude and the lag of the thermal tide. 

This problem should be conceptually analogous to the classical problem of a periodically heated 
slab (3). With that in mind we derive the thermal response of an idealized slab and will later show that 
the frequency dependence of the solution captures, rather surprisingly, striking features of the thermal 
tide response. 

Let us consider a slab of heat capacity C per unit area at temperature T, radiating like a black body, 
and receiving an energy flux F =Fo + 5F. Let us assume that the temperature response can be linearized 
around a mean temperature, T = To + 5T. The temperature evolution can thus be written 

^ Cj^ST = Fo-OsBT^ + 8F-4osBTi6T. ( 6 ) 


Now, if we assume that the mean temperature is such that the slab is in mean thermal equilibrium 
(Fq = OsbTq^, and that the perturbation is periodic of the form (because we consider only one 
component of the Fourier decomposition of the insolation, e.g. the semi-diurnal one), the complex 
amplitude of the response is given by 


ioCdT = dF -4osBTidT ^ dT 


5F 1 
2coqC l + io/{2coo) ’ 


(7) 


with (Do = 2 o^bTo /C, which is twice the inverse of the usual thermal timescale. Interestingly, if the 
heat capacity of the atmosphere/surface system to be modeled can be estimated (for example by C = 
CpPs/g + Qurf). one can have an order of magnitude estimate for coo. 


2 < 756 ^ 0 ^ 

CpPs/,? T Csurf 


( 8 ) 


Note that Tq is the emission, and not surface, temperature. This order of magnitude estimate is actually 
confirmed by the numerical model. If further scalings are needed, let us note that (Oq'xTq oc ^3/4 p^j. 
a given star we can further say that ft)o oc n. 

By analogy with our problem where the important frequency is a = 2ft) — 2n, let us model the 
pressure bulge amplitude and lag by 


— 


_ ^0 _ 

1 +i{(0 — n)/(Oo 


Im(^a) = qo 


{(0-n)/(0Q 
\ + [{(0-n)/(0Qf' 


(9) 


where qo{pfi,Fo) and 0 )o{ps,Fo) can be fitted from a set of numerical atmosphere models. Despite 
its simplicity, and the lack of proper dynamical treatment. Fig. 1 shows that this parametrization of 
the torque actually captures salient features of the complex frequency dependence of the atmospheric 
response, especially when considering the huge uncertainties in other parameters. Then, thanks to the 
simplicity of our approach, we can describe the thermal response of a given atmosphere with only two 
parameters that can be computed for a given set of atmospheric parameters (mean pressure, incoming 
flux, etc.). Results for our various sets of simulations encompassing the habitable zone can be found in 
Table 1. 

Giving an analytical order of magnitude estimate of the amplitude of the disturbance (especially 
in terms of pressure) is however rather difficult because one would need to solve both dynamical and 
radiative transfer equations. The important, and counter-intuitive, differences between our results for 
the Earth and Venus show that more work is needed to find a reliable analytical estimate. For the time 
being, the model thus still needs to be calibrated on numerical results as the one presented in Table 1. 

The main difference with earlier analytical treatments (4, 5) is that we account for the cooling feed¬ 
back of the atmosphere (last term in Eq. (2)) which yields a response oc 1/(mg + /a) instead of l/ia. 
The former formula was found in ref. 3, but their equivalent of cOq was unknown and assumed to be zero, 
yielding the later formula. Our treatment thus avoids any discontinuity in the zero frequency limit and 
self-consistently predicts a linearly increasing, then saturating, lag angle. 
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5 Existence and calculation of asynchronous equilibrium spin states for 
various rheological models 

As shown in ref. 5, independently of the tidal models used, finding equilibrium spin states reduces to 
solving the following equation 


f{(0-n) 


b^{2(0 — 2n) 
bg{2(0 — 2n) 


K, 


( 10 ) 


Here we take the convention that the obliquity is always 0° and that ft) is negative for a retrograde 
rotation, but we keep in mind that for each equilibrium state (ti),0°) found, the (—ft), 180°) state is also 
an equilibrium (5). 

The location and number of solutions, however, does depend on the specific model used. If fhe 
resfricfion / of / fo fhe positive frequency range is monofonic around possible equilibria, only one syn¬ 
chronous and fwo asynchronous spin slates exisl (5). An imporlanl qualilafive novelty of fhe presenl 
work is lhal ba is nol monofonic. As a resull, more equilibria polenlially exisl, depending on fhe rhe¬ 
ology used lo describe fhe bulk of fhe planel. In general, fhis rheology is paramefrized by a frequency 
dependenf quality factor Q{o) wifh 


f7g(a) = sign(a) ^||^|j , (11) 

where k 2 is fhe Love number of degree 2. 

Recenl sludies suggesl lhal lidal dissipation inside ferreslrial planels can be accuralely described by 
fhe so called Andrade rheology which accounls for bolh fhe viscous behavior of fhe manlle in fhe zero 
frequency limil (very close lo synchronization) and fhe low frequency dependence af higher frequencies 
(similar lo fhe consfanf-Q model). Thus, we firsl discuss fhe exislence of equilibrium asynchronous 
spin sfales wifh fhe Andrade model. Because a simple analylical solufion is difficull lo attain, we Ihen 
show how fhe qualilafive behavior of fhe solufion can be underslood in terms of fwo simpler models (fhe 
conslanl time lag and conslanl-Q models) lhal allow us to find analylical solulions in limil cases. 

We do nol consider Ihe viscoelastic Maxwell model (which disregards inelastic behaviors which are 
imporlanl in Ihe body of Ihe planel) because ils frequency dependence al higher frequencies 
1 /a) is to steep to property explain existing dala for lerreslrial planels (28). Indeed, Ihe Earlh’s Q 
changes by slighlly more lhan an order of magnilude belween Ihe Chandler wobble period (aboul 440 
days) and seismic periods of a few seconds (28, 30, 31) whereas Ihe Maxwell model would predicl a 
change of aboul 7 orders of magnilude ( 1 440days/1 sec | ss 4 x 10^). 

5.1 The Andrade model 

Allhough viscosity dominates Ihe dissipation in Ihe manlle al very small frequencies (a <C where 
Tm is the Maxwell time), inelasticity takes over at higher frequencies. Both these behaviors are encom¬ 
passed in the Andrade model (see ref. 28 for details) whose dissipation is described by 

kijo) ^ _ g sin(g^) r(l + «)) __ 

2(^) 2 ('('I +A2)a-f gi^“%“cos(^)r(l -f a))^-f -f gi-“T^“sin(^)r(l -ha))^ 

where A 2 is a constant of order unity depending on the mass, radius and rigidity of the planet, and 
r is the gamma function. For the Earth, the Maxwell time, Tm, is about 500 yr, and the index of the 
frequency power law at higher frequencies, a, is estimated between 0.15 and 0.3 (we will use 0.2 as a 
fiducial value). 

As is visible in Fig. 2, when Ihe amplilude of Ihe Ihermal tide reaches a Ihreshold, Ihere are five 
differenl equilibria. This is due to Ihe facl lhal / is non monofonic over Ihe range of frequencies of 
inleresl. However, to give an analytical formula for Ihis Ihreshold we need to simplify Ihe problem 
furlher. 
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5.2 The constant-2 model 


An important simplification comes from the fact that the Maxwell time for terrestrial planet is much 
larger than the forcing periods we are envisioning. Therefore, for forcing periods of the order of days to 
months (a » the tidal dissipation scales as a^“ which is weakly frequency dependent (a ~ 0.2). 
This is needed to explain the low frequency dependence of the Earth Q (see above) which stems from 
the low frequency dependence of inelastic processes in rocky mantles. Originally, it is why it has very 
often been assumed that 2(a) is independent of frequency for rocky planets, giving rise to the so-called 
constant-2 (or constant phase lag) model. 

Indeed, as seen in Fig. 2, this very simple model reproduces most of the features of the full Andrade 
model. In particular, we recover the existence of four asynchronous equilibria when a given criterion is 
met. This criterion is given by 




qoQ 



> 1 . 


( 12 ) 


We can also calculate the location of the various spin states 


(o^=n±a)Q(jQ + e^Yl-l^, (13) 

where e = 1 for stable equilibria and e = — 1 for unstable ones. This equation can be deduced from Eq. 
7 in ref. 5 applied to our specific rheology. It is noteworthy that, because of the shape of the various 
torques, the degree of asynchronism {\( 0 ± — n\) does not continuously goes to zero (as can be seen in 
Fig. S3), and is always greater than coq. 

Because it is completely analytical, and considering the small differences with the full Andrade 
model compared to the uncertainties in the various parameters involved, we mainly discuss the results 
of this model in the body of the article. 

The main limitation of this model is that it predicts a discontinuity near synchronization where the 
frequency vanishes. It thus cannot be used to describe what happens to the synchronous spin state. 
Various authors have tried to overcome this problem by changing Q near a = 0 so as to have a linear 
behavior but the range of frequency over which they made the transition was mostly unconstrained. This 
is resolved in the Andrade model. 


5.3 The viscous model 


Indeed, in the zero frequency limit, the Andrade model tells us that the mantle acts as a viscous fluid, 
for which the time delay between the perturbation and the response is constant so that Q ^ ^ At G. In 

this regime, / = k 2 Atwo i+{(j/(oof ' specific model, / being monotonic, there is at most 

two asynchronous spin states which exist only when the following criterion is met 


7vis — 


5 K, 


qo 


2471 Ko k2At(Oo 


> 1 . 


(14) 


Then, the two equilibrium rotation rates are given by = n ± G)o\/7vis ~ 1- It can be shown that this 
criterion is also sufficient for the synchronous spin state to be an unstable equilibrium. 

A simple estimation for Venus yields /vis < 10^^. This is because for a rocky, Earth-like mantle, 
Tm is extremely large (~ 500 yr), and k2At = 3/2 x A 2 TM in this context. As a result, for a body whose 
rheology obeys the Andrade law, it is very likely that the synchronous spin state is a stable equilibrium. 
This has important implications as discussed below. 

For sake of completeness, note that, if the planet were to behave as a viscous fluid over the whole 
frequency domain with a much lower viscosity, then the situation might be a little different. Calibrating 
the time lag on the Earth Moon system, one finds At ~ 630^ <C Xm (29). Then, criterion (14) would be 
much easier to fulfill, and the planet could have two stable, asynchronous equilibria and one unstable 
synchronous one, as depicted in Fig. S4. Although often used for its linear behavior, data showing the 
inadequacy of this model to describe solid planets are accumulating (2S). We thus mention this only to 
help the reader make the link with the various models available in the literature. 
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5.4 Critical distance for the existence of asynchronous equilibrium spin states 

Rewriting the criterion given by Eq. (12), one can see that planets are expected to have non-synchronous, 
equilibrium spin states when located beyond a critical orbital distance given by 





7 \ 1/3 



(^gRpM. 



1/3 


(15) 


which is plotted in Fig. 3 for an Earth-like planet with Ps = 1 and 10 bar and the value of <70 corre¬ 
sponding to F = 1366W.m^^. The only dependence in the flux received by the planet appears through 
qo. However, although qo increases with flux as expected, it changes by less than ~ 25% across the 
habitable zone despite a flux variation of a factor 3 (see Table 1). As a result, as shown in Fig. S5, for a 
given pressure, the explicit dependence of the critical distance on the flux received by the planet is small 
compared to the pressure dependence (and the uncertainty on the bodily tides). For planets outside the 
habitable zone, receiving much larger/lower fluxes, additional simulations would be needed. 


6 Stellar properties and habitable zone 

As discussed in the main text, we often have to relate the flux received by a given planet, its orbital 
distance and the stellar mass. For example, in Fig. S3, once the atmospheric properties and the substellar 
top-of-atmosphere insolation F are chosen, the orbital distance will vary with the stellar mass. To do so, 
the bolometric luminosity (L*) of a star of M* is given by (32) 




(16) 


where jx = log(M*/M 0 ). Then the distance at which a planet receives F is given by a = y 

For the boundaries of the habitable zone shown in Fig. 3 and used throughout the text, we used the 
parametrization of the outer edge from ref. 24. The parametrization of the inner edge is taken from the 
same reference but is rescaled to take into account recent results from 3D climate models (20); i.e. the 
dependence on the effective temperature of the star is left unchanged, but the inner edge is set at 0.95 AU 
around a Sun-like star. 
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7 Trapping in synchronous resonance and constraints on atmospheric 
evolution 

If a solid planet with a non-zero triaxiality ever approaches near the 1:1 synchronous spin-orbit reso¬ 
nance, the conservative restoration torque exerted by the star on its permanent asymmetry can, in some 
cases, trap the planet in spin-orbit resonance. As pointed out in ref. 5 and 30, even if the synchronous 
rotation rate were made unstable by thermal tides, spin-orbit trapping could, in some cases, stabilize a 
planet in this state if it ever went near it (if the atmosphere was formed a long time after the planet’s 
formation, for example). 

However, our more realistic treatment of both thermal and body tides has made this issue rather 
moot. Indeed, with a realistic rheology, the synchronous spin state is found to be stable (see Fig. 2) so 
that a planet whose spin approaches this region will end-up being synchronized even without spin-orbit 
trapping. 

But we know that Venus did not get caught in a synchronous spin state while having an internal 
structure quite similar to the Earth (which is well modeled by an Andrade rheology). So this shows that 
planets can, and some will, avoid evolving toward such a stable synchronous state. 

This creates a very interesting, new constraint. With our treatment, even without resonance trap¬ 
ping, if Venus spins down fast enough to have a rotation rate between the two unstable equilibria in 
Fig. 2 (which could act as repelling barriers) before the atmosphere is created, then even once the at¬ 
mosphere is outgassed, the spin-down will probably proceed toward synchronization. This tells us that 
the atmospheric build-up on Venus had to occur on a shorter timescale than the tidal despinning, which 
constrains both the dissipation efficiency in the planet and its atmospheric evolution. 

Fortunately, this is in good agreement with the recent developments in the theory of atmosphere for¬ 
mation in general, and on Venus in particular. These developments suggest that just after the formation 
of the planet itself, the mantle is in a magma ocean state and the atmosphere consist of tens of bars of 
CO 2 and water in solution equilibrium with the magma. Then the magma ocean solidifies (affer one 
fo fens of Myr) and hundreds of bars of CO 2 (up fo fhousands for H 2 O) are released fo fhe afmosphere 
because fhese gases are poorly soluble in solid rocks (33). 

Affer fhaf, for Venus, ref. 34 suggesfs fhaf fhe afmosphere never cooled down sufficienlly fo allow 
fhe wafer fo condense and fhaf fhe wafer jusf escaped fo space in fhe firsf fens fo hundreds Myr living 
a fhick CO 2 afmosphere behind. Thus, even if fhis afmosphere subsequenfly evolved, if is unlikely fhaf 
if ever gof orders of magnifude fhiner fhan fhe presenf one. The afmospheric fide may fhus have been 
efficienf during mosf of Venus’ history. 

For fhe Earfh, wafer is expecfed fo have condensed when fhe surface cooled enough, creating fhe 
oceans and reducing fhe mass of fhe afmosphere fo a value fhaf is commensurate wifh fhe currenf one. 
Then, fhe longer, secondary, afmospheric evolufion began. Wifh fhaf in mind, if is unlikely fhaf fhe Earfh 
ever had a much fhinner afmosphere fhan ifs presenf one. Some dafa even indicate fhaf if could have been 
fhicker during some periods in ifs evolufion (35). Again, fhis supporfs fhe idea fhaf fhermal fides may 
have been efficienf for mosf of fhe hisfory of planefs in fhe habifable zone who could refain a significanf 
afmosphere. 
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Figure SI: Surface pressure anomaly created by the thermal tide. A-D: Snapshots of the spatial 
distribution of the departure of the surface pressure from its mean value created by the thermal tide. 
E-H: Spatial distribution of the semi-diurnal component only. The location and direction of motion of 
the substellar point is shown with a white arrow. From top to bottom, the length of the solar day is 
decreased from 200 (planet near synchronization) to 25 days for a model with a 10 bar mean surface 
pressure on a 100 d orbit. This shows the increase in the lag and decrease in amplitude with increasing 
forcing frequency. 
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Figure S2: Amplitude and phase of the thermal tide for the Earth model. A: Map of the amplitude 
of the semi-diurnal component of the surface pressure tide (normalized to its maximum value) for the 
Earth validation model. B: Amplitude of the diurnal tide normalized to the maximum amplitude of the 
semi-diurnal component. This shows that, on Earth, the semi-diurnal component is dominating. C: 
Phase of the semi-diurnal tide. The average value at the equator is ss 150°, close to the expected value. 
To be compared with Fig. 2S.3 of ref. 1. 
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Figure S3: Degree of asynchronicity as a function of stellar mass for Earth-like planets in the 
habitable zone. |co± — u|/u is computed from Eq. (13) for planets with different atmospheric masses 
(characterized by the surface pressure, p^), near both the inner and the outer edge of the habitable zone. 
We only show stable equilibria. 
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Figure S4: Equilibrium spin states of the planet in the weak friction approximation. Atmospheric 
(dashed), gravitational (dotted) and total (solid) torque as a function of spin rate for the viscous model. 
Arrows show the sense of spin evolution. A: weak atmospheric torque, only one equilibrium, syn¬ 
chronous spin state exists (blue disk). B: bifurcation point (a = ac). C: the atmospheric torque is strong 
enough to generate two stable, asynchronous, equilibrium spin states (blue disks; one is retrograde in 
the case shown). The synchronous spin state is unstable. 
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Figure S5: Dependence of critical "non-synchronous" orbital distance on planetary flux. For each 
pressure {p^ = 1 and 10bar), the critical distance, ac, beyond which planets are expected to exhibit 
non-synchronous equilibrium spin states is shown for two extreme values of qo representative of the 
conditions in the inner (red) and outer (blue) regions of the habitable zone. The dashed line corresponds 
to the pure CO 2 outer edge case. Other features are defined in Fig. 3. For a given atmospheric pressure, 
the critical distance varies little with the value of the flux received by the planet in the atmospheric 
model, especially when compared to the variations with pressure. 
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